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' We study the time evolution of the reduced density matrix of a system of spin- 1/2 particles 

interacting with an environment of spin-1/2 particles. The initial state of the composite 
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system is taken to be a product state of a pure state of the system and a pure state of the 
environment. The latter pure state is prepared such that it represents the environment at a 
given finite temperature in the canonical ensemble. The state of the composite system evolves 

\q ' according to the time-dependent Schrodinger equation, the interaction creating entanglement 

between the system and the environment. It is shown that independent of the strength of 

£NJ ' the interaction and the initial temperature of the environment, all the eigenvalues of the 

reduced density matrix converge to their stationary values, implying that also the entropy 
■ of the system relaxes to a stationary value. We demonstrate that the difference between the 

canonical density matrix and the reduced density matrix in the stationary state increases 
as the initial temperature of the environment decreases. As our numerical simulations are 
' necessarily restricted to a modest number of spin-1/2 particles (< 36), but do not rely on 

d ' time-averaging of observables nor on the assumption that the coupling between system and 

environment is weak, they suggest that the stationary state of the system directly follows 
from the time evolution of a pure state of the composite system, even if the size of the latter 
cannot be regarded as being close to the thermodynamic limit. 
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Statistical mechanics is one of the cornerstones of modern physics 1, 2 ) but its foundations 
are still subject of much research. 3 ~ 23 ) Some fundamental questions, such as how the canonical 
distribution emerges from the interaction between a system S and its environment E, have only 
been partially resolved. Answers to this question are important since the canonical ensemble 
is widely used to calculate the thermodynamic quantities of a system at a given temperature. 

As well for classical as quantum systems it is well-known that if the interaction between 
a system S and its much larger environment E having a large number of degrees of freedom 
and a dense distribution of energy levels, is weak, the system S is described by a canonical 
ensemble when the composite system S + E is described by the microcanonical ensemble with 
a given total energy. Usually, the derivation of the canonical distribution is discussed under 
the hypothesis that each state in the microcanonical ensemble has equal probability. 1 ' 2 ) In the 
case of quantum systems, it has recently been shown that the microcanonical mixed state for 
the composite system S + E is not a required starting point for the system S to be described 
by a canonical ensemble, but that S + E being initially in a randomly picked pure state with 
small energy fluctuations is sufficient. 7,10 ' 11,14 ) The characteristic that even if the state of the 
composite quantum system S + E corresponds to a single wave function only, the reduced 
density matrix of S is canonical for the overwhelming majority of wave functions in the 
subspace corresponding to the energy interval encompassed by the microcanonical ensemble, 
is referred to as canonical typicality after Ref. 10 ) Not explicitly mentioning canonical typicality, 
this characteristic had already been used to calculate the density of states (DOS) of quantum 
many-body systems. 24,25 ** More recently it has been shown that canonical typicality has a 
classical counterpart: For typical probability distributions defined on an energy shell of the 
classical composite system S+E, i.e. not necessarily microcanonical distributions, the marginal 
probability distribution corresponding to the system S exhibits the canonical form. 26 ) 

In this paper we focus on the equilibration obtained from the dynamics of relatively small 
closed quantum systems (containing less than 36 spin-1/2 particles) that we compose from a 
small system S and a much larger but still relatively small (containing less than 32 spin-1/2 
particles) environment E. We use general quantum spin-1/2 Hamiltonians to describe S and 
E and we do not put any restriction on their energy spectra. We study the conditions under 
which the stationary state of the system S is represented by a canonical ensemble density 
matrix, although the standard conditions, such as the environment E being very large and 
the coupling between S and E being weak, are not (necessarily) fulfilled. The approach we 
use is to first solve the time-dependent Schrodinger equation (TDSE) of the composite system 
S + E of spin-1/2 particles numerically and then analyze the behavior of the reduced density 
matrix of the system S, obtained by tracing out the environment E, which was initially 
prepared in a randomly picked "typical" pure state with a given temperature. 
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Earlier work adopting this approach showed that a nanoscale environment prepared in 
a uniform random superposition of all states (corresponding to an environment at infinite 
temperature) can drive the system to the state with a canonical distribution 27 ) and elucidated 
the effect of frustration and connectivity on decoherence and relaxation processes. 28_3 °) In this 
paper, we extend the approach to study the decoherence and relaxation properties of nanoscale 
magnets embedded in a nanoscale magnetic environment at finite temperature. 

Our simulation results show that, independent of the strength of the interaction between S 
and E and the initial temperature of E, S evolves to a stationary state of which the properties 
strongly depend on the initial temperature of E. This equilibration is remarkable given the 
relative small size of E, since usually in equilibration studies the hypothesis of having a 
large environment is essential. 14, 20 ) We show that for sufficiently large initial temperatures of 
E, the stationary state of S is represented by a canonical ensemble density matrix at some 
finite effective temperature. For decreasing temperatures, the reduced density matrix of S 
deviates from the canonical density matrix. The deviation increases for decreasing values of 
the interaction strength between S and E. 

The paper is organized as follows. In Section 2, we discuss the model, define the quantities 
of interest and summarize the essentials of the simulation method used. The simulation results 
are presented in Section 3. A discussion and conclusion is given in Section 4. 

2. Generalities 



In general, the state of a closed quantum system is described by a density matrix. 31, 32 ) 
The canonical ensemble is characterized by a density matrix that is diagonal with respect 
to the eigenstates of the Hamiltonian H, the diagonal elements taking the form exp(— (3Ei) 
where j3 = l/ksT is proportional to the inverse temperature (kg is Boltzmann's constant and 
is taken to be one in this paper) and the E^s denote the eigenenergies of H. 1 ' 2 ^ 

The time evolution of a closed quantum system is governed by the TDSE. 31,32 ) If the 
initial density matrix of an isolated quantum system is non-diagonal then, according to the 
time evolution dictated by the TDSE, it remains non-diagonal and the quantum system never 
approaches the thermal equilibrium state with the canonical distribution. Therefore, in order 
to equilibrate the system S, it is necessary to have the system S interact with an environment 
E, also called a heat bath. Thus, the Hamiltonian of the composite system S + E takes the 
form H = H$ + He + Hse, where H$ and He are the system and environment Hamiltonian, 
respectively and H$e describes the interaction between the system and environment. 
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2.1 Model 

To study the evolution to the canonical ensemble state in detail, we consider a general 
quantum spin- 1/2 model defined by the Hamiltonian H = Hg + He + H$e where 

ns—X ns 

Hs = - E E J i,i S i S ^ C 1 ) 

i=l j=i+la=x.y,z 
n E -I n E 

H *> = E E E W?. ( 2 ) 

i=l j=i+l a=x,j,2 

^ = -EE E a&w- ( 3 ) 

i=l 7=1 <x=x,y,z 

Here S 1 and / denote the spin-1/2 operators of the system and environment, respectively 
(we use units such that H and are one). The total number of spins in the system and 
environment are denoted by ns and tie, respectively. 

The spins of the system are arranged in a ring and interact via a isotropic Heisenberg 
interaction Jf- = J. The spins of the environment are all connected with each other and with 
all the spins of the system. Previous work 27-30 ) has shown that it is expedient, though not 
essential to take for the spin-spin interactions A"-, and Qfj uniform random numbers in the 
range [— |A| , |A|], and [— |f2| , respectively. Relative to other choices of these interactions, 
the randomness of the interaction parameters A and Q and the high connectivity of the spins 
generally reduce the decoherence and relaxation time to reach the stationary state of the 
reduced density matrix. Note that we do not put any restriction on the energy spectra of the 
Hamiltonians describing S and E. 

2.2 Initial state 

We prepare the state of the system S and of the environment E separately at t < and 
then bring them in contact with each other at t = 0. Specifically, we construct the initial pure 
state of the composite system S + E, p(Qi) = \^(0))(^(0)\ where 

'* (0)) " |S) (« B |.-'".|«')V (4) 

with \$e) = Ci\4>i) denoting the state of the environment with the coefficients q generated 
randomly according to the prescription given in Ref. 25 ) and {|</>i)} being an orthonormal set 
of basis states which, in our simulation software, are the usual direct products of the spin up 
and down states. Numerically, the imaginary-time propagation by e~^ HE ^ 2 is performed by 
means of a Chebyshev polynomial algorithm. 33 ~ 36 ) To prepare the environment in its ground 
state (/3 = oo), we use the standard Lanczos method. 37 ) 

It follows directly from Ref. 25 ^ that for any observable = 0) of the environment 

(*(0)\X E (t = 0)|*(0)) w Trp E X E (t = 0), (5) 
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the approximation improving as the inverse square root of the dimension of the Hilbert space 
of the environment (see Appendix). Therefore, we may consider the state |^(0)) as "typical" 
in the sense that if we measure observables of the environment, their expectation values 
agree with those obtained from the canonical distribution of the environment at the inverse 
temperature (3. Note that in practice, it is often sufficient to consider only one random state 
\&e) (see Appendix). 

The assumption of random phases in the initial state has been instrumental in the deriva- 
tion of the quantum master equation, 38, 39 ) a key equation in the theory of non-equilibrium 
statistical mechanics. Within the quantum master equation approach, the approach to equi- 
lbrium of a quantum system is well understood. 1, 39 ** Although there may be an apparent 
similarity with the use of the random initial states that we use in the present work, there 
is no relation between the random initial states and the random phase assumption in the 
derivation of the master equation. In the present work, random initial states are a convenient 
computational device only: As we show in the Appendix, their use effectively eliminates the 
need to compute traces of operators and allows us to work with pure states only. Below, we 
also demonstrate explicitly that the use of random initial states is not essential for the main 
conclusions of this paper by starting the simulation from the initial state with all spins up. 

2.3 Time evolution 

A pure state of the composite system S + E evolves in time according to (in units of h = 1) 

D a D E 

|tf(t)> = e- UH \*(0)) = ^2 y £c( l ,p,t)\i,p), (6) 

i=l p=X 

where the states are just another notation of the complete set of orthonormal states 

in the spin-up - spin-down basis and D$ = 2™ 3 and De = 2™ E denote the dimension of the 
Hilbert space of the system and environment, respectively. 

Numerically, the real-time propagation by e~ ttH is carried out by means of the Chebyshev 
polynomial algorithm, 33 ~ 36 ) thereby solving the TDSE for the composite system starting from 
the initial state |^(0)). This algorithm yields results that are very accurate (close to machine 
precision), independent of the time step used. 40 ) 

2.4 Reduced density matrix 

The state of the quantum system S is described by the reduced density matrix 

p(t) = Tr E p(t), (7) 

where p (t) is the density matrix of the composite system at time t and Tr^ denotes the trace 
over the degrees of freedom of the environment. The system S is in the canonical state if the 
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reduced density matrix takes the form 

m = er fiHs /Tr s e-^ t (8) 

where Trg denotes the trace over the degrees of freedom of the system S. In terms of the 
expansion coefficients c(i,p,t), the matrix element of the reduced density matrix reads 

D E D E 

p id {t) = Tr E^2^c*(i,q,t)c(j,p,t)\j,p)(i,q\ 

p=l q=l 



De 



^c*(i,p,t)c(j,p,t). 

p=l 



(9) 



2.5 Data analysis 



We analyze the time-dependent data of the reduced density matrix in various ways. First, 
at each time step (in units of r = tt/W), we diagonalize the (non-negative definite) reduced 
density matrix itself and study the time-dependence of its eigenvalues. We define the variance 
of the set of eigenvalues at t and tf by 



var(t) = 



5>i(t) - \i(t f )y 

i=l 



(10) 



where Xi(t) is the ith eigenvalue of p(t). Usually, tf is taken to be the final time of the 
simulation. From the eigenvalues, we also compute the entropy of the system 

D S 

S(t) = -Tvp{t)\np{t) = -^Aj(t)ln \(t). 

i=l 

We characterize the degree of decoherence of the system by 

a(t) 



(11) 



D S -1 D S 

E E \p^)\\ (12) 
\ i=i j=i+i 

where Pij(t) is the matrix element of the reduced density matrix p in the representation 
that diagonalizes Hg. Clearly, a{t) is a global measure for the size of the off-diagonal terms of 
the reduced density matrix in the representation that diagonalizes Hg- If c{t) = the system 
is in a state of full decoherence (relative to the representation that diagonalizes H$)- 

Assuming that the system S, evolving in time according to the TDSE, relaxes to the 
canonical state we expect that p(t) ~ p(b) for t > to where to is some finite time and b 
denotes the effective inverse temperature of S. The difference between the state p(t) and the 
canonical distribution p(b(t)) is conveniently characterized by 



5(t) 



Ds 



Ds 



^E(fti(i)-e- 6 W Ei /E rKf)Ei 



(13) 



i=l 
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Fig. 1. Simulation results for var(i/r) and S(t/r) (see Eqs. (10) and (11)) and <x(t/r) and S(t/r) (see 
Eqs. (12) and (13)) for different values of the initial temperature of the environment, as obtained 
by solving the TDSE for a system consisting of a ring of four spins coupled to an environment of 
eighteen spins. The spins of the environment are all connected with each other and with the four 
spins of the system. The initial state of the system S is given by | titi) an d the initial state of the 
environment E is given by a random pure state. Model parameters: J = —0.5, A = 0.2, CI = 0.5, 
tf/r = 190; Sample interval r = 7r/10. Solid line: j3 = 1; Dotted line: /3 = 3; Dashed line: j3 = 6; 
Dash-dotted line: (3 = oo. 



with 

b{t) = = . (14) 

If the system relaxes to its canonical distribution both 5(t) and a(t) are expected to vanish, 
b(t) converging to the effective inverse temperature b. 

For any function /(.) of the system Hamiltonian H$, we define the averages 

(f(H S )) m = Trp(t)f(Hg), (15) 

and 

(f(H s )) b = Tre- bH sf(H s )/Tre- bH s. (16) 
Then, application of the Schwarz inequality yields 

\(f(H s )) m -(f(Hs)) b \ 2 < 5 2 (t)Trf(H s ), (17) 
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Table I. Data for er, S, Sp and Ep, taken at the last time step of the simulation run. Sb and Eb are 
calculated according to Eq. (16). 



p 


b 


a 


5 


Sb 


Sp 


E b 


Ep 






J = 


-0.5, A 


= 0.2, n 


= 0.5 






1 


0.807 


0.003 


0.003 


2.704 


2.706 


-0.165 


-0.162 


3 


1.851 


0.011 


0.015 


2.392 


2.412 


-0.400 


-0.388 


6 


2.526 


0.019 


0.027 


2.080 


2.114 


-0.543 


-0.526 


oo 


3.044 


0.031 


0.030 


1.816 


1.852 


-0.638 


-0.621 



showing that the deviations of the energy and entropy of the system from their values in the 
canonical ensemble vanish linearly or faster with 5(t). 

3. Results 

Most of our simulations have been carried out for systems consisting of four spins coupled 
to an environment of eighteen spins. We have verified that our conclusions do not depend on 
details such as the connectivity of the spins in the environment or the size of the composite 
system by simulating triangular lattices, regular square lattices and so on with up to 35 spins 
(data not shown) . 

In Fig. 1, we present the simulation results for var(i/r), S(t/r) (see Eqs. (10) and (11)), 
a(t/r) and S(t/r) (see Eqs. (12) and (13)) for the case that there is a fairly strong coupling 
between the system S and the environment E (|A/J| = 0.4) and for different values of 
the initial temperature of the environment. From Fig. 1, it is clear that, independent of the 
initial temperature of the environment, all the eigenvalues Aj(i) of the reduced density matrix 
converge to a stationary value. This implies that also the entropy of the system S approaches a 
stationary value, suggesting that the system S relaxes to an equilibrium state. We emphasize 
that the data presented in this paper are obtained without any averaging procedure other 
than the one intrinsic to quantum theory. 

From Fig. 1, it follows from the data for a(t) that if j3 = 1, the reduced density matrix 
of the system S converges to a diagonal matrix in the representation that diagonalizes H$, in 
other words, the system S has lost virtually all coherence {a — > 0) and as 5 — > with time, 
the system S relaxes to the canonical system. The same holds for (3 = 0, see Ref. 27 ) With 
increasing /3, the difference between the reduced density matrix and the canonical distribution 
(for the system defined by H$) increases slightly, as indicated by the fact that the values of 
a and 5 increase with /3. 

It is instructive to analyze more quantitatively, the data taken at the end of these par- 
ticular simulation runs. In Table I we collect the results for the various quantities of interest. 
As f3 increases, a increases too, indicating that the deviation from a diagonal matrix (with 



8/21 



J. Phys. Soc. Jpn. 



Full Paper 



0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 




1 

\ 
1 





























50 



100 

t/T 



150 



200 




100 

t/T 



150 



200 



3 
2.5 
2 

S. 1.5 

CO 

1 

0.5 





50 



100 
t/T 



150 



200 



50 



100 

t/T 



150 



200 



Fig. 2. Simulation results for var(i/r), S(t/r), a(t/r), and 6(t/r) for different values of the energy 
scale of the couplings between the spins in the environment. The data are obtained by solving 
the TDSE for a system consisting of a ring of four spins coupled to an environment of eighteen 
spins. The spins of the environment are all connected with each other and with the four spins 
of the system. The initial state of the system S is given by | tltl) arL d the initial state of the 
environment E is given by a random pure state. Model parameters: J = —0.5, A = 0.2, /? = 6, 
tf/r= 190; Sample interval r = 7r/10. Solid line: ft = 0.125; Dotted line: f2 = 0.25; Long-dashed 
line: fl = 0.5; Dash-dotted line: f2 = 0.75. Short-dashed line: fi = 1. 



respect to the basis that diagonalizes Hs) increases, merely reflecting the fact that the de- 
coherence processes become less effective as the temperature decreases. Nevertheless, even at 
zero temperature, the difference between the reduced density matrix and the canonical density 
matrix is quite small, of the order af a few percent, and so are the differences for the entropy 
and energy. Thus, it seems that even for fairly strong coupling between the system and the 
environment (|A/J| = 0.4), the nanoscale environment drives the even smaller system to a 
state that, within a few percent, is described by the canonical state of the system, albeit with 
an effective temperature that does not agree with that of the environment (compare f3 and b 
in Table I). 

3. 1 Energy scale of the environment 

To investigate the effect of the parameter O, the energy scale of the states of the environ- 
ment, we performed simulations for = 0.125,0.25,0.5,0.75, 1. In Fig. 2 we present data for 
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Fig. 3. Same as Fig. 1 except that the system-environment interaction strength A = 0.02. Model 
parameters: J = -0.5, A = 0.02, fl = 0.5, tf/r = 11000; Sample interval r = tt/10. Solid line: 
(3=1; Dashed line: (3 = 6. 



(3 = 6. From earlier wor 

k 28-30) for p = and p 

= co we know that the more the ranee and 
the width of the energy spectrum of the environment and the system match each other, the 
better the system decoheres. Taking this into account we can easily understand the behavior 
of a: When the value of O decreases, the value of a decreases. 

Except for £1 < \ J\ (recall J = —0.5 in this paper), both a(t) and 5{t) relax to fairy small 
values, the difference between the reduced density matrix at t/r = 200 with the canonical 
states being of the order of a few percent (data not shown). For Q > |J|, the qualitative 
behavior is different: Although S(t), a(t) and 5(t) to relax to their stationary values, the dif- 
ference between the reduced density matrix and the canonical state is significant, as indicated 
by 8{t/r = 200) ~ 0.1. Qualitatively, this can be understood as follows. Keeping the number 
of spins of the environment and the range of S — E interactions A fixed, increasing f2 increases 
the spectral range of the environment, that is the spacing between the energy levels of the 
environment increases. This effectively reduces the mixing of the eigenstates of the system and 
the environment by the S — E interactions, which in turn leads to a reduction of the effect of 
decoherence by the environment and the probabilities for the system to exchange energy with 
the environment. 
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Table II. Same as Table I except that the coupling A = 0.02 instead of A 



Full Paper 
= 0.2. 



f3 b a 5 Sb Sp Et> Ep 





J = - 


-0.5, A 


= 0.02, 


n = 0.5 






1 0.874 
6 3.349 


0.003 
0.047 


0.005 
0.207 


2.692 
1.659 


2.696 
2.125 


-0.180 
-0.687 


-0.175 
-0.495 



3.2 Weak interaction between system and environment 

For the simulation data presented earlier, the system-environment interaction strength A 
was chosen such that |A/J| = |A/f2| = Oil) which is far away from the weak coupling regime. 
If we reduce A, we may expect that the time scale of decoherence and relaxation processes 
increases with A. To keep the amount of computer time required to reach the stationary state 
within reasonable limits, we have chosen to reduce A by a factor ten, that is we take A = 0.02 
and consider this to be the "weak coupling" case. The simulation results for this parameter 
choice are presented in Fig. 3. Note that compared to Fig. 1 and Fig. 2, the time scale to 
reach the stationary state has increased by a factor of about one hundred. For /3 = 1, there is 
no qualitative change compared to the case of A = 0.2: The state of the system converges to 
the canonical state. 

However, for (3 = 6 the relatively large values of a(t) and 5(t) signal that the decoher- 
ence process is not very effective and that the deviation from the canonical distribution is 
significant, as is shown quantitatively in Table II. Nevertheless, the eigenvalues of the reduced 
density matrix relax to their stationary values and so does the system entropy. Qualitatively, 
this can be understood by the same argument as the one used at the end of Section 3.1. 
Keeping the number of spins of the environment and the spectral range of the environment 
fixed, decreasing the range of S — E interactions A effectively reduces the mixing of the 
eigenstates of the system and the environment, which in turn leads to a reduction of the effect 
of decoherence by the environment and the probabilities for the system to exchange energy 
with the environment. 

3.3 Properties of the stationary state 

The observation that the eigenvalues of the reduced density matrix, and therefore also 
the entropy of the system, approach stationary values for sufficiently long times seems to be 
generic, independent of the values of the parameters A/| J| < 1, Q/\ J\ < 1, (3 and the initial 
state of the system S or the environment E, as illustrated in Fig. 4. As a matter of fact, in 
our large collection of simulation results (most results not shown) there is no evidence of the 
contrary. However, the stationary state itself depends on all the above mentioned parameters. 
The time scale on which the equilibration occurs strongly depends on the system-environment 
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Fig. 4. The entropy of the system as a function of t/r for the same system as in Fig. 1 except 
that the initial state of the system and environment is the state with all spins up (solid line) or 
| *(())) = e _,3ifB / 2 |$)/($|e~^ ffE |$) 1 / 2 (dashed line) where |$) is a uniform random superposition 
of all states of the whole system and /3 = 6. 
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Fig. 5. Left:The difference dis(t/r, Ai/r) = Tr\p(t) — p(t + At)\/2 between the reduced density ma- 
trices p(t) and pit + A)t as a function of t/r for a fixed value At/r = 1000, extracted from the 
simulations that yield the data presented in Fig. 3. Solid line: /3 = 1; Dashed line: (3 = 6. Right: 
The energy of the system as a function of t/r for a time interval chosen in which there is exchange 
of energy between the system and the environment, showing that the dynamical evolution of the 
system is non-Markovian. 



interaction strength A. As expected, reducing A increases the time scale of decoherence and 
relaxation processes. Although our simulation results are obtained for a very small quantum 
system connected to a relatively small environment, they are in good agreement with findings 
for interacting large quantum systems that seem to evolve in such a way that any small 
subsystem equilibrates, under the condition that the Hamiltonian has no degenerate energy 
gaps and that the state of the composite system contains sufficiently many eigenstates. 20 ) 

However, our simulation results also show that the eigenstates of the reduced density 
matrix generally do not evolve to a stationary state. In Fig. 5 (left) we show representative 
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results of a measure for the distance, the trace distance dis(i/r, Ai/r) = Tr\p(t) — p(t + Ai|/2, 
between the two reduced density matrices pit) and p(t + At) as a function of t for a fixed value 
of At, as obtained from the simulations that yield the data presented in Fig. 3. For /3 = 1, 
the case in which the difference between the reduced density matrix and the canonical state 
is small, dis(t/r, Ai/r) becomes very small for t/r > 8000, suggesting that the eigenvectors of 
the reduced density matrix also converge to their stationary values (in concert with the fact 
that in this case, the reduced density matrix is very close to the canonical state). In contrast, 
for /3 = 6 and t/r large, dis(f/r, Ai/r) levels off at a value which is not small. In other words, 
although the eigenvalues of the reduced density matrix relax to their stationary values, the 
eigenvectors of the reduced density matrix exhibit nontrivial quantum dynamics, even though 
the entropy of the system has reached a stationary, nonzero, value. 

In the case of a single spin and two-spins interacting with an environment, the dynamics 
of eigenvectors of the reduced density matrix in the stationary-state was observed earlier 
through quantum oscillations of spin expectation values. 41, 42 ) In Fig. 5(right), we present 
results for the energy of the subsystem, as obtained from the simulation that produced the 
data of Fig. 3(right). It is clear that the system energy mostly decreases, with energy going 
from the system to the environment. However, for some time intervals (one being shown), 
the system energy increase, indicating that the environment transfers energy to the system, 
a characteristic feature of non-Markovian processes. 43 ' 44 ) 

4. Conclusion 

We have presented a simulation study of a small magnetic system coupled to a nanoscale 
magnetic environment. The quantum dynamical evolution of the composite system was ob- 
tained from the direct numerical solution of the TDSE. Our analysis, albeit numerical, does 
not involve any approximation, does not rely on time- averaging of observables nor does it 
assume that the coupling between system and environment is weak. 

The most striking result of our analysis is that all the eigenvalues of the reduced density 
matrix relax to stationary values, implying that the entropy of the system relaxes to a station- 
ary value. In this sense, the nanoscale environment drives the system to a thermodynamically 
stationary state, a feature which is usually attributed to macroscopic environments only. 

Furthermore, we have shown that under suitable but fairly general conditions, the reduced 
density matrix in the stationary state is close to the canonical state of the system, albeit not 
with the same temperature as the one of the environment. The difference between these two 
states generally increases as the temperature of the environment decreases. 

For /3 ~ 0, the initial state of the environment has all the features of a "canonical typical- 
ity" state and qualitatively, our results are in concert with the theory of canonical typicality: 
The reduced density matrix relaxes to the canonical state. 10 ' 11 ) As (3 increases the concept 
of "canonical typicality" no longer applies in the strict sense, as explained in the Appendix. 
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Yet, qualitatively, we may interpret our findings in terms of the effective dimension d e J^ of 
the environment: 11 ) If /3 ~ 0, <fj = O(De), which is a large number for the systems that we 
have used in our simulations. As (3 increases, the number of states of the environment effec- 
tively available for decoherence and relaxation decreases (see also the results of the density of 
states in Fig. A-l). As this implies that d e J^ decreases, it is to be expected that the difference 
between the reduced density matrix and the canonical state increases. n ) 
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Appendix 

The use of random initial states has played a central role in developing "fast" (i.e. 0{D)) 
algorithms to compute the density of states and other similar quantities. An early application 
of such an algorithm to electron motion in disordered alloy models was given by Alben et 
al. 45 ) It was shown that the eigenvalue spectrum of a particle moving in continuum space can 
be computed in the same manner. 46 ) Fast algorithms of this kind proved useful to study vari- 
ous aspects of localization of waves, 47-49 ) other one-particle problems 35, 50 > 51 ) and many-body 
problems. 24,25 ) The rigorous proof that this approach has remarkable statistical properties, 
namely that for large D the statistical error vanishes as 1/VT), was given in Ref. 25 ) 

Following Ref., 25 ) we consider real random variables x\, y\, . . ., xd, yD, taking values in 
the interval [—00, +00] and distributed according to the probability density 

/(xi,yi,.. .,x D ,y D ) = -^TT 5 ( x i + ■■■ + VD ~ !); ( A>1 ) 

where T(D) is the Gamma function. Writing c n = x n + iy n for n = 1, . . . , D, we construct the 
random state 

D 

|$)=^c„|^), (A-2) 

n=l 

where is a complete set of orthonormal basis states of the D-dimensional Hilbert space, 

which for the derivation that follows need not be specified further. 
From Eq. (A-l), it directly follows that 

/+00 
c k f(x 1 ,...,y D )dx 1 ...dy D = 0, (A-3) 
-00 



{c* k c k ,) = I c* k c k tf(xu...,y D )dx 1 ...dy D = S k ,k' D » ( A-4 ) 
J —00 

and that (c k Ck>) = 0. It also follows from Eq. (A-l) that !<!>) is a unit random vector with a 
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uniform probability density on the hypersphere of dimension D — 1. 
Next we consider the projected state 

D 

m/3/2)) = e- m/2 \^>) = dje-^ /2 \Ej), (A-5) 

i=i 

where Ej (\Ej)) denotes the j-th eigenvalue (eigenstate) of the Hamiltonian H and 

D 

dj = ^RjnCn, (A-6) 

i=i 

where Rj n = (Ej\ip n ) is the unitary transformation matrix. The probability density of the 
random variables dj reads 

fid!, ...,d D ) = ^^(|di| 2 + • • • + \d D \ 2 - 1). (A-7) 

Hence the dj are distributed uniformly over the D-dimensional hypersphere. Furthermore, we 
have (dj) = and (d*jdji\ = D^ x 8jji. 
Normalizing the state Eq. (A-5) yields 

D , —/3E/2 d 



Zf=iWe~^ i=i 



where 



and 



aj = - j d £L , (A-9) 
E,=i \ d j\ 2 Pj 



Pi = _ flB! , , (A-10) 



is the Boltzmann weight for the state j. In general, the probability density of the coefficients 



dj is not uniform. 



Next we want to show that for sufficiently large D, we may replace EjLi \dj\ 2 Pj by its 
average over the distribution Eq. (A-7), that is by D^ 1 . To prove this, we compute the average 
of 

X 2 = (d- 1 -J2\dj\ 2 p)j , (A-ll) 
with respect to the distribution Eq. (A-7). We have 

(x 2 ) = ^ 2 -2^x:^<l^l 2 ) + EEw<KfK'l 2 > 

j=y i=i i'=i 

D D 

i=i j'=i 
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Using Eq. (2), (A12) and (A23) in Ref., 25 ) we find 

/\rl l 2 U l 2 \ - 25 33' | 1 ~ S Jf 1 + S Jf (ft -iq\ 

W \*s>\ ) - d(dTT) + mp + i) ~ D(D + iy (A ' 13) 



yielding 



[X 2\ _ V-„2 



D(D + l)^ D*(D + 1) 

D - 1 1 , . 

Invoking Markov's inequality, 52 ) it follows that 

P(X 2 > < D" 1 . (A-15) 

In words, the probability that the error X 2 is larger than D -1 is smaller than D~ l . For the 
case at hand D increases exponentially with the number of spins. For instance for a system 
of 18 spins and a random state with probability density Eq. (A-l) 

P(X 2 > 0.38 x 10~ 5 ) < 0.38 x 10 -5 , (A- 16) 

suggesting that for all practical purposes, it is safe to assume that X ~ and that with 
probability very close to one, the projected state Eq. (A-8) can be written as 

D 

\$(P/2))=D 1 / 2 J2d jP 1 / 2 \E j ). (A-17) 



1/2 

Putting bj = djPj , the probability density of random variables bj is 



m,..,M=^( M +---+ M -i)^(n^) (A-i8) 

From these calculations, the following conclusions can be drawn about the projected state 
Eq. (A-17): 

(1) From Eq. (A-17) and the properties of random variables {dj}, it follows directly that on 
average 

w/m\HPm) = (a-iq) 

assuming that D is sufficiently large. 25 ) Note that if (<S>\Ei) ^ (\E\) denoting the non- 
degenerate eigenstate) we have lim^ tX) ($( ( S/2)|F|$(/3/2)) = (Ei\Y\Ei), independent of 
D. 

(2) According to Eq. (A- 18), in general the projected state Eq. (A-17) is randomly distributed 
on a hyper-ellipsoid, not on a hyper-sphere, in the Hilbert space. This suggests that it may 
be of interest to extend the concept of canonical typicality from states on an hyper-sphere 
to states on an hyper-ellipsoid by introducing a non-zero inverse temperature beta. 

For completeness, we briefly discuss simulation results for the case that we prepare the 
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Tabic A-l. Same as Table I except that the composite system S + E is prepared in a random state 
that is typical for the composite system being in the canonical state at inverse temperature (3. 






b 


a 


5 


s b 


Sfi 


Sp 


E b 




Ep 








J = 


-0.5, A 


= 0.2, 


n = 0.5 








1 


0.976 


0.004 


0.002 


2.671 


2.666 


2.671 


-0.202 


-0.208 


-0.202 


3 


2.602 


0.025 


0.016 


2.042 


1.838 


2.003 


-0.557 


-0.630 


-0.569 


6 


3.744 


0.045 


0.077 


1.462 


0.633 


1.210 


-0.742 


-0.920 


-0.799 


oo 


4.575 


0.029 


0.132 


1.092 


0.000 


0.673 


-0.832 


-1.000 


-0.908 



Table A-2. Same as Table A-l except that the coupling A = 0.02 instead of A = 0.2. 



p 


b 


a 


6 


s b 




Sp 


E b 


Ep 


Ep 








J=- 


-0.5, A 


= 0.02, 


n = 0.5 








1 


0.998 


0.003 


0.001 


2.666 


2.666 


2.666 


-0.207 


-0.208 


-0.207 


6 


5.981 


0.032 


0.020 


0.637 


0.633 


0.573 


-0.920 


-0.920 


-0.928 



25 




Fig. A-l. Simulation results for the density of states (DOS) of the environment, conditional on the 
initial state of the environment, for the model parameters of Fig. 3 corresponding to the weak- 
coupling case. Solid line: j3 = 1; Dotted line: /3 = 3; Dashed line: j3 = 6. 



composite system S + E, not just the environment E, in a pure state that is typical for the 
canonical state at a given /3, namely 

|lp(0)> = <5FH5)W (A ' 20 > 

where H is the Hamiltonian of the composite system, |3>) = ^ ■ Cj\cpj) is a pure state of the 
composite system, Cj is generated randomly according to the prescription given in Ref., 25 ) and 
{|</>j)} is an orthonormal set of basis states. According to the theoretical analysis presented 
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earlier in this Appendix, averages taken with this pure state will yield the values that are 
equal to those taken with respect to the canonical state with temperature T = 1//3. 

Solving the TDSE with the initial state Eq. (A-20) yields results for cr(t), S(t), var(i) 
and S(t) that, up to small fluctuations, are constant in time, apparently consistent with the 



taken as an indication that the state Eq. (A-20) is also "typical" in the sense of "canonical 
typicality". However, for the same reasons as those given earlier, this conclusion would be 
incorrect unless (3 — > 0. 

Obviously, the nanoscale models that we study are not "thermodynamic" in the usual 
sense, but the number of states is quite large (Fig. A-l shows the DOS of ~ 4 x 10 6 states) 
and it is possible, in principle, to find many (but not macroscopically many) states with 
energies in a narrow interval such that the DOS in this interval is almost constant. However, 
this is not sufficient to apply conventional statistical mechanics arguments: What is required 
is that also the number of particles (22-35 in our work) is large. Otherwise the fluctuations of 
the energy (which are proportional to the inverse square root of the number of particles) are 
not small and the equivalence between microcanonical and canonical ensemble is no longer 
guaranteed. 1 ' 2 ) 

Our simulations show that the dependence of <j(t), 5(t), var(t) and S(t) on j3 is qualitatively 
the same as in the case that we use the product state Eq. (4) as the initial state: In all cases 
considered, the eigenvalues of the reduced density matrix converge to stationary values. For 
comparison with Table I, in Tables A-l and A-2 we give the data extracted from the simulations 
using Eq. (A-20) as the initial state. It is clear that both Tables show the same qualitative 
features as a function of f3. For f3 = 1 and A = 0.02 our results are in concert with Tasaki's 
analytical results 7 -* for weak coupling (/3A <C 1 in the notation of Ref. 7 )). 

When we prepare the composite system in a pure state that is typical for the canonical 
state at a given /3, the simulation data of cr(t), S(t), var(t) and S(t) show very little time- 
dependence (as discussed above). However, as suggested by the data in Tables A-l and A-2, 
the reduced density matrix is not equal to e~@ Hs /Trse~^ Hs but is "renormalized" by the 
interaction H$e as can be seen by the perturbative treatment that follows. 

Up to the second order in the interaction Hamiltonian H$e we have 53 - 1 



statement in Ref., 10 ) that u p 



,*(*) 



pa even at t = for typical wave functions" . This may be 



e -P(Hs+H E +H SE ) 



e 



P{H S +H E ) 



dx e -(P-x)(Hs+H E ) HsEe -x(H s+ H E ) 



+ I dx ^ dy e ~{P--)(Hs+H E ) HsEe -(x-y)(H s +H E ) HsEe --y(H s+ H B ) 



Jo J 




18/21 



J. Phys. Soc. Jpn. Full Paper 

yielding Z = Z S Z E {1 -zi + z 2 ) where Z s = Tr s e~^ Hs , Z E = Tr E e~^ HE , 

— — Tr f dxe-^-^ Hs+H ^H SE e-^ Hs+H ^ 
sZ E J 



ZsZ E 

P 

Z,sZ E 



Tre-P {Hs+HE) H SE , (A-22) 



and 

z 2 = —!—Tr f dx T ( iye-^- a; )^ + ^)^5Ee- (:!; - s/){Hs+ ^ ) ^£;e- y ^ s+HE) , (A-23) 
ZsZ E J Jo 

are the first and second-order correction, respectively. Up to second order in the interaction 
Hamiltonian H$ E , the reduced density matrix therefore reads 



m 



Tr s Tr E e-P( H s+ H B+H SE ) ~ Z S Z E (1 - Zl + z 2 ) 

/ dxTr E e-P HE e xHs H SE e 
Jo 



p\ X + Zx -z 2 -\z\- I dx Tr E e-P HE e xH *Hs E e- x 



Z E 

+ J dxj*dy Tr E e xHs e-^- x+y)HE H SE e- {x - y){Hs+HE) H SE e- yHs ^j 

+0(H 3 SE ), (A-24) 

showing that in the thermal equilibrium state, due to the interaction, we should expect a 
deviation from the canonical distribution of the system. 

A direct numerical calculation of the various contributions is beyond our current capabili- 
ties and we therefore leave this calculation for future research. However, comparing the energy 
and entropy of the system in the canonical state with the corresponding values obtained from 
the simulation of the composite system provides some idea of how much the reduced density 
matrix changes as a result of the interaction H$ E - From Table A-l, we conclude that for 
A = 0.2 the differences \Sp — Sp\ and \Ep — Ep\ are of the order of 10% or more, except 
for (3 = 1. Hence A = 0.2 definitely does not correspond to the case of weak interaction. 
From Table A-2, it is clear that A = 0.02 and /3 = 1 correspond to weak interaction between 
environment and system because Sp w Sp and Ep Ep but for /3 = 6, \Sp — Sp\ is of the order 
of 10%, hence not small. Thus, we conclude that for both A = 0.02 and A = 0.2, the effect 
of the interaction on the reduced density matrix is significant if /3 = 6, even if we prepare the 
composite system in a typical canonical state. 
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